Screening properties of Gaussian electrolyte models, 
with application to dissipative particle dynamics 



Patrick B. Warren, 1 ']^] Andrey Vlasov, 2 Lucian Anton, 3 >n and Andrew J. Masters 3 

1 Unilever R&D Port Sunlight, Quarry Road East, Bebington, Wirral, CH63 3JW, UK. 
2 Department of Chemistry, St. Petersburg State University, 
26 Universitetsky prosp., 198504 St. Petersburg, Russia. 
3 School of Chemical Engineering and Analytical Science, 
University of Manchester, Manchester, M13 9PL, UK. 
(Dated: March 4, 2013) 

We investigate the screening properties of Gaussian charge models of electrolyte solutions by 
analysing the asymptotic behaviour of the pair distribution functions. We use a combination of 
Monte-Carlo simulations with the hyper-netted chain integral equation closure, and the random 
phase approximation, to establish the conditions under which a screening length is well defined and 
the extent to which it matches the expected Debye length. For practical applications, for example in 
dissipative particle dynamics, we are able to summarise our results in succinct rules-of-thumb which 
can be used for mesoscale modeling of electrolyte solutions. We thereby establish a solid foundation 
for future work, such as the systematic incorporation of specific ion effects. 



I. INTRODUCTION 

Dissipative particle dynamics (DPD) has seen 
widespread uptake in modelling soft condensed matter 
[U [2]. The attractions are obvious: by coarse grain- 
ing over the atomistic degrees of freedom one can ac- 
cess the relevant length and time scales with only mod- 
est computing requirements. Polymer phase behaviour 
[3 H] j polymer dynamics [5 , polymer rheology [6j , sur- 
factant mesophase formation kinetics [7j, the properties 
of amphiphilic bilayers [HI IH| , and the properties of col- 
loidal suspensions [10] have all been investigated by the 
method. 

Charged systems such as anionic and cationic sur- 
factants, water-soluble polyelectrolytes, charge-stabilised 
colloidal suspensions, and mixtures of these [TT] , form a 
large subclass of widespread practical importance. In 
these systems there is often the requirement to model 
the supporting electrolyte. This can be done implic- 
itly, for example with the Poisson-Boltzmann equation, 
or explictly by incorporating ions as charged particles 
in the simulation. In the latter case, particularly for 
DPD where soft interactions are the norm, it is natural 
to smear the point charges into charge clouds. The di- 
vergence of the long-range Coulomb law as r — > (where 
r is the center-center separation) is replaced by a smooth 
cutoff, thus ensuring thermodynamic stability according 
to a theorem by Fisher and Ruelle [T2] . 

The precise form of the charge smearing is often tuned 
to the choice of numerical algorithm and a consensus on 
the best approach has yet to emerge. Groot introduced 
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a grid-based method with linear charge smearing [11) . 
Later Gonzalez-Mclchor et al. examined an Ewald-based 
method with exponential charge smearing [13] . Here 
we study a related Ewald method with Gaussian charge 
smearing. This choice can be used to simplify the Ewald 
algorithm, and connects with recent work on the so-called 
ultrasoft restricted primitive model (URPM) [Mlfl6] . In 
principle the differences between smearing methods can 
be subsumed into short-range part of the interparticle 
potential, though the details are the subject of ongoing 
investigations. 



To study the screening properties of our Gaussian elec- 
trolyte model, we use a combination of Monte-Carlo 
(MC) simulations, the hyper-netted chain (HNC) integral 
equation closure, and the random phase approximation 
(RPA), to analyse the asymptotic behaviour of the pair 
distribution functions. The programme is as follows. In 
the next two sections we define the mesoscale electrolyte 
model and the tools used to analyse it. We then present 
results demonstrating that, for typical applications, HNC 
can be relied upon to deliver accurate results (it is no 
exaggeration to say that HNC is up to ten million times 
faster than MC). We then use HNC to explore the screen- 
ing properties of the model, establishing the conditions 
under which a screening length is well defined (i. e. on the 
low density side of a Kirkwood line in the phase diagram) 
and the extent to which the screening length matches 
the expected Debye length. We further establish the do- 
main of applicability of the much simpler RPA, which 
gives relatively simple expressions for the Kirkwood line 
and the screening length. We emphasise that our ap- 
proach could easily be applied to other smeared charge 
electrolyte models. Mindful of this, and the utility of a 
fast, accurate, multicomponent HNC solver in general, 
we have released our FORTRAN 90 HNC code as fully 
documented open source software |17j . 



2 



II. MODEL 

We now describe the Gaussian charge model for elec- 
trolyte solutions. The potential energy is given by a sum 
of pairwise terms, split into short range and long range 
(electrostatic) contributions, 

' X ''<••• U H= U ii + U ii- (!) 

The short range piece is given by 

/3U*= { ^(i-^iAc) 2 (r«<r c ) 
I (nj > r c ) 

and the long range piece is given by 

^6-^ -*(£)■ < 3 » 

In these /3 = l/k^T is the inverse of the temperature T 
measured in units of Boltzmann's constant fcs, T%j is the 
centre-centre separation between particles i and j, 
is a dimensionless short range repulsion amplitude which 
depends on the particle types, Ib is the Bjerrum length 
which plays the role of an electrostatic coupling constant, 
Zi and Zj are the valencies measured in units of an elemen- 
tary charge, and r c and a are length scales which mea- 
sure, respectively, the range of short range repulsion and 
the size of the Gaussian charge cloud. The short range 
part of the potential corresponds to the standard DPD 
interaction law [3J. The long range part corresponds to 
the interaction between Gaussian smeared charges with a 
radial charge distribution (27rcr 2 ) -3 / 2 exp(— r 2 /2a 2 ). The 
function erf(r/2<r) — > rjia^pK) as r — > 0, thus ensuring 
the Coulombic divergence is replaced by a smooth cutoff. 

We will consider up to three species of particles, cor- 
responding to positively and negatively charged ions of 
valencies z + and z_ at densities p+ and p_ , and a third 
neutral solvent species at a density po. The total ion den- 
sity will be denoted by p z = p+ + p_ . The total overall 
density will be denoted by p = p + p z . In the case where 
there is no solvent, p = and p = p z . We adopt the 
convention that the valency includes the sign as well as 
the magnitude. Overall charge neutrality then requires 
z + p + + z_p_ = 0. We do not necessarily suppose the va- 
lencies are of the same magnitude. We shall label species 
by Greek indices, a, f3 — (0, +, — ). The system volume is 
V and thermal averages will be denoted by (•). 

We first consider a special case. The aforementioned 
URPM is an unsolvated equimolar mixture of Gaussian 
charge clouds, corresponding to the choice Ay = 0, 
po = 0, and \z±\ = 1. The URPM is governed by a 
dimensionless density, p z o 3 , and a dimensionless cou- 
pling constant, Ib/o~, which plays the role of an inverse 
temperature. The model exhibits marked clustering for 
Ib/o- ^ 30, and a condensation transition for Zb/ct > 100, 
for densities in the range p z a 3 0.01-0.1 (these esti- 
mates are translated from the results shown in Fig. 5 



in Ref. 14). The physics behind the phase transition re- 
mains somewhat unclear [16(118). but the phenomenology 
can be viewed as a reflection of stability issue for point 
charges mentioned in the introduction. It quantifies the 
onset of the 'danger zone' as the point charge limit is 
approached. To avoid these artefacts the implication is 
that we should attempt to keep Ib/o~ < 30. However 
for practical applications there is already a strong incen- 
tive to make a as large as possible, to reduce the cost of 
computing the electrostatic interactions. Usually this is 
enough to ensure that low temperature URPM artefacts 
are avoided. 

In the general case the properties of the model are gov- 
erned by three length scales, r c , cr and Ib, the repulsion 
amplitude matrix Aij , the choice of valencies z± , and the 
densities p z and p. The parameter space is thus poten- 
tially very large. Our strategy to reduce the complexity 
is to consider the mapping to the underlying atomistic 
system. This requires us to distinguish between physi- 
cal units in which the length scales and densities are ex- 
pressed in SI units; and simulation units in which length 
scales and densities are expressed in units of r c or a. 

In standard DPD the choice pr 3 = 3 is usually made 
|3J, and we will adopt the same here. In addition one 
usually introduces the notion of a 'mapping number' 
N m , giving the number of solvent molecules represented 
by one DPD solvent particle. Given this, the value 
of r c in physical units is determined by the identity 
pN m V m / Na = 1, where V m is the solvent molar volume 
and Na is Avogadro's number [T3]. If water is the solvent 
(V m — 18 x 10~ 6 molm -3 ), and with the conventional 
choice N m — 3, one has in physical units r c = 0.645 nm. 

Next consider the Bjerrum length. In physical units 
this is Ib — e 2 /(47re r eofcB7 1 ) where e is the elementary 
charge, e r is the relative permittivity of the solvent, and 
e is the permittivity of free space. For water at room 
temperature Ib ~ 0.7 nm ss 1.09 r c . Since a z:z elec- 
trolyte is equivalent to a 1:1 electrolyte with Ib increased 
by a factor z 2 , there is considerable interest in exploring 
higher values of Ib- In the present work we shall explore 
up to Zb = 10 r c ( « 7nm in physical units), which covers 
many cases of interest. 

For an electrolyte at a molar concentration c s , the mi- 
croscopic ion density is 10 3 c s Na- This is readily con- 
verted to a simulation density by multiplying by r 3 .. For 
example a 0.1 M 1:1 electrolyte solution would be repre- 
sented by p z r 3 = 0.032 (note that p z counts both species 
of ion) . To cover the typical range of electrolyte concen- 
trations we therefore consider p z r 3 in the range 10 _3 -1. 

The above considerations do not yet impinge on the 
choice of a. This is a central theme of the present study 
and will be discussed extensively below. 

Finally we discuss the repulsion amplitude matrix. In 
the present study we will only consider a constant repul- 
sion amplitude matrix Ay = A, leaving the extension to 
unequal repulsion amplitudes for future work. We use 
either A = 25 motivated by standard DPD [3J, or A = 
corresponding to the situation in the absence of short 
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range repulsions. In the latter case, of course, it does not 
make sense to include the neutral solvent species since 
it would just form an ideal gas in the background. It 
has been suggested that A should be chosen to match 
the solvent compressibility However this introduces 
the danger if A is too large one will encounter an order- 
disorder transition driven by the short range repulsions. 
A less rigorous criterion is to demand only that the sol- 
vent be relatively incompressible, so that d(J3p) /dp ^> 1 
where p is the pressure. This is satisfied by A = 25, for 
which the solvent is clearly still a liquid with only moder- 
ate structure. Moreover much work has been done based 
on this value, which we will therefore continue to use. 

To summarise, in the remainder of this work we shall 
consider mainly two classes of models: either the pure 
URPM comprising unsolvated Gaussian charges, or the 
'solvated' case containing in addition a neutral species 
and short range repulsions between all particles. (It is 
possible to consider an intermediate case where short 
range repulsions are added to the URPM, however this 
does not generate any new insights.) As already men- 
tioned the URPM is characterised by the dimensionless 
density p z a 3 (with \z±\ = 1 implying p± = p z /2) and 
coupling strength Ib/c. The length scale r c plays no role, 
except, perhaps, as a 'fiducial' length. On the other hand 
the solvated case (with the 'standard' choice pr 3 = 3 and 
A = 25) is characterised by the dimensionless densities 
P± r c: an d dimensionless ratios Ib/i" c and cr/r c , where all 
except o~/r c are fixed by the mapping to the underlying 
physical system. 

From a practical point of view a ^ r c is commonplace, 
and it shall be important to pay attention to the units 
of length when mapping between the solvated case and 
the pure URPM. In the text we shall endeavour to be 
always explicit about this, and where the figure annota- 
tions use implicit units we shall always state the choice 
of units in the caption. A symmetric z:z electrolyte in 
the solvated case can be mapped to the URPM with a 
renormalised Ib — > z 2 Ib- An asymmetric electrolyte can- 
not be mapped onto the URPM but we shall discuss this 
case only rather briefly. It is worth bearing in mind that 
it is quite straightforward to apply the tools developed 
here to all these cases. 

III. TOOLS 

A. Pair distribution functions and screening 

Given the model is governed solely by pair interac- 
tions, the thermodynamic properties are completely de- 
termined by the pair distribution functions g a p (r) . In ad- 
dition the screening properties are also determined by the 
asymptotic behaviour of these functions. Specifically, the 
screening length A features in the asymptotic behaviour 
of the total correlation functions, 

e -r/X 

h a p{r) = g a p(r) - 1 ~ (r ->• oo) (4) 



provided the asymptotic decay is purely exponential. The 
decay length defined in this way is unique to each state 
point and does not depend on the identity of the species 
of charged particles under consideration. 

As the density decreases the screening length ap- 
proaches the Debye-Hiickel limiting law behaviour, A — > 
Ad, where the Debye length is 

Ad = («bE«4)" 1/2 . (5) 

Conversely, as the density increases the screening length 
gets smaller but there comes a point where the asymp- 
totic decay of the total correlation functions ceases to be 
purely exponential and instead becomes damped oscilla- 
tory. This transition defines a line in the phase diagram 
known in charged systems as the Kirkwood line [15], or 
more generally a Fisher- Widom line [217] . 

For applications, one would hope that the actual 
screening length will hew as closely as possible to the 
expected Debye length (at least, as long as the latter 
is well defined). The extent to which this can be made 
so is the central theme of the present work. For exam- 
ple, a 1:1 electrolyte at 0.1 M concentration has a Debye 
length Ad ~ 0.96 nm rj 1.5 r c . If we simulate this in the 
present model with the choice a = r c , the asymptotic 
decay of the total correlation functions would be oscilla- 
tory and we would not even be on the right side of the 
Kirkwood line. On the other hand if we use a = r c /2 
the asymptotic decay would be purely exponential with 
A/Ad = 0.94, thus the actual screening length would be 
only 6% different from the Debye length. This example 
will be worked through in more detail below. 

B. Integral equation theory 

Given that the underlying electrolyte model is a fluid 
mixture, it is natural to think of using multicomponent 
integral equation theory to calculate the structural and 
thermodynamic properties 2 1 2~i1. Further, since the 
interactions are soft, one expects that the hyper- netted 
chain (HNC) integral equation closure should work well. 
We find this is indeed the case. We also find that for 
parameters typical of 1:1 electrolytes, the random phase 
approximation (RPA) also works well. 

The starting point is the multicomponent Ornstein- 
Zernike (OZ) relation which defines the direct correlation 
functions c a p(r). In reciprocal space the OZ relation is 

h a /3 = Ca/3 + J2f Pi Caj h"yf3 (6) 

where the spatial Fourier transform of a function x(r) is 
defined by x(k) = Jd 3 r e~ lkr x(r). The HNC closure is 
defined in real space, and is 

h a p = exp(-/3C/ Q/3 + h a p - c a0 ) - 1 (7) 

where U a p is the pair potential between particles of 
species a and /3. The solution of these coupled equa- 
tions is numerically quite demanding and in the present 
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case we exploit the accelerated convergence schemes orig- 
inally proposed by Ng [22H24] . We typically solve the 
distribution functions on a grid of size 4096 points at a 
grid spacing 0.01 r c , so that the functions are calculated 
out to r « 40 r c where all trace of structure has typically 
vanished below the numerical precision of the calculation. 
We find, however, that the schemes fail to converge for 
Ib/p > 10. This loss of solution has also been observed 
by Coslovich, Hansen and Kahl [T5], and may be indica- 
tive of a mathematical property of the HNC rather than 
a numerical problem. 

Numerically much less demanding is the RPA closure, 
which is given by 



C a /3 



-pu a 







(8) 



Because there are no hard cores, the RPA is the same 
as the mean spherical approximation (MSA). Unlike the 
HNC, the RPA can be solved for all values of Ib/ct 
although it may yield unphysical results (for example 

hap < -1). 

The pressure p and the internal energy density (U) /V, 
can be solved from the pair functions [STJ [23] . The pres- 
sure can be found either by the virial or compressibility 
routes. These do not give exactly the same result because 
the HNC closure breaks thermodynamic consistency. In 
practice, for the present applications, we have found the 
two routes differ by typically at most a few percent. Spe- 
cific results for the RPA thermodynamics can be found 
in Refs.HHandlH] 



C. RPA solution of the URPM 

Coslovich, Hansen and Kahl solve the RPA for the 
URPM [15] . and we have recently revisited the prob- 
lem in terms of the low temperature phase behaviour 
[18] . The relevant properties of the RPA solution are de- 
scribed here. URPM symmetry implies that in the RPA 
the correlation functions are given by h±±(r) — ±h(r). 
Inserting this into the OZ equations, with the RPA clo- 
sure, reveals 



h(k) = 



~4ttI b exp(—a 2 k 2 ) 
k 2 + k 2 ^ cxp(-er 2 fc 2 ) 



(9) 



where k 2 ^ = AttIbPz is the square of the Debye wavevector 
(i. e. fc D = 1/A D ). 

The real space total correlation functions can (in 
principle) be obtained by expressing the Fourier back- 
transform of Eq. ([9]) as a contour integral in the complex 
/c-plane. The behaviour of the correlation functions is 
therefore determined by the poles of h{k) in the upper 
half plane. As a particular consequence, the asymptotic 
behaviour of h(r) as r — > 00 is determined by the po- 
sition^) of the pole(s) closest to the real axis [To] 125], 
There are two cases. If the nearest pole to the real axis 
is purely imaginary, the asymptotic behaviour of h(r) is 



purely exponential, with a decay length set by the dis- 
tance of the pole from the real axis. Alternatively if the 
nearest poles to the real axis are complex, the asymp- 
totic behaviour is damped oscillatory. Clearly, then, the 
Kirkwood line is determined by the crossover between 
these two scenarios, typically when a purely imaginary 
pole nearest to the real axis collides with the next-nearest 
pole, to form a complex pair which subsequently move off 
the imaginary axis. 

In the present case, writing q — ak and = ct&d, the 
poles of h(q) are determined by the solutions of 



Z 2 exp(q 2 ) 



-A-kIbPzO- 



(10) 



These solutions can be expressed in terms of the Lambert 
W function which solves We w = z. For the asymptotic 
behaviour of h(r) the most relevant solution is given by 
q 2 = Wo(— ?d) where Wo is the principal branch of the 
Lambert W function [26l|27]. If q£ < 1/e, W (-qk) is 
negative real and the corresponding poles of Eq. ([9]) (in 
the complex q-plane) are at q = ±z|H / o(— q 2 ^)] 1 ^ 2 . The 
corresponding decay length is given by 



Arpa = ax \W a (-4irl B Pz(T 2 



-1/2 



(11) 



Note that Arpa — > Aq (from below) as p z — > 0. If q 2 ^ > 
1/e, Wo(— 3d) is complex and the asymptotic decay of 
h(r) is damped oscillatory. We therefore identify q 2 ^ 



1/e as the Kirkwood line [16]. Equivalently, Eq. (11) 
requires 



AnelBPzV 2 < 1 , 



(12) 



with equality determining the location of the RPA Kirk- 
wood line. For general applications the charge density 
in Eqs. (11) and (12) should be taken to be the ionic 
strength, defined by p z — J2 a z aPa- The RPA can of 
course be solved for the solvated URPM case. We leave 
discussion of this to a separate publication. 



D. Monte-Carlo methods 

We benchmark HNC against MC simulations. We use 
an NVT ensemble with standard single particle trial dis- 
placements in the usual Metropolis scheme [T]. We cal- 
culate the energy of a configuration from U ~ U s + U L , 
where II s = J2 l>3 &om Eq. (§, and U L = £ t>J U% 
from Eq. The latter is re-expressed as an Ewald sum, 



27T/ E 



]>> fc |Qi< 



k<k c 



2(7 



(13) 



where Ak = k~ 2 exp(— (J 2 k 2 ), and Qk = 
the reciprocal space charge density. This is just the stan- 
dard Ewald result omitting the real space contribution 
[TJ [28] . The last term is a self energy correction; this can 
of course be omitted from the MC acceptance criterion, 
but it is essential to retain this correction when making 
comparisons with integral equation theory. 
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FIG. 1. (color online) Pair distribution functions for the 
URPM at two state points on opposite sides of the Kirkwood 
line (see text for details), plotted as |r/i±±| versus r to illus- 
trate the asymptotic behaviour. Lines are HNC, data points 
with error bars are MC. Lengths are expressed in units of a. 



The first term in Eq. ( 13 ) is a sum over a discrete set 



of wavevectors, commensurate with the simulation box 
dimensions, such that k = |k| < k c where the cut-off 
is chosen so that exp(— cr 2 k 2 ) is sufficiently small. For 
the simulations reported below we use ak c = 4. There 
are about 47rfc 3 /3 x (2w) 3 /V discrete wavevectors in the 
sum (the second factor is the density of wavevectors in 
reciprocal space). This means that the computational 
cost of evaluating the sum varies as 1/cr 3 . As a practical 
consideration, this is a strong motivation for making a 
as large as possible. 

The only other point to make about the Ewald imple- 
mentation concerns the calculation of the pressure, 



l>j 



+ 



2nl B 
3V 2 



Y / Ak\Q k \ 2 (l~2a 2 k 2 ) 



(14) 



k<k c 



The first and second terms in this are the ideal gas re- 
sult and the standard virial result for pair interactions. 
The third term follows from Eq. ( 13 ) by calculating 
-(dU h /dV) [25]. Obviously the form of this term means 
that it can be easily evaluated alongside the energy. 



IV. RESULTS 

A. Comparison between HNC and MC 

We first consider the URPM as a baseline. Fig. [I] 
shows the pair distribution functions plotted as \r h(r)\ at 
two state points on either side of the Kirkwood line: (a) 
p z a 3 = 0.02 and Ib/& = 1 where the asymptotic decay 
is purely exponential; and (b) p z a 3 — 0.2 and Ib/ct = 10 
where the decay is damped oscillatory. The agreement 
between HNC and MC is excellent. These MC simula- 
tions are expensive to perform even in the absence of a 
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FIG. 2. (color online) URPM thermodynamics along two 
isotherms showing — p cx (solid HNC lines with square MC 
data points) and — {U)/SV (dashed HNC lines with diamond 
MC data points) as a function of density. The Debye-Hiickel 
limiting law is shown as a dotted line in the two cases. Lengths 
and densities are expressed in units of a, and thermodynamic 
quantities in units of k^T/a 3 . 
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FIG. 3. (color online) Pair distribution functions for a sol- 
vated model at the indicated state point. Lines are HNC, 
data points with error bars are MC. From top to bottom the 
curves are: — ; go+ = go- ~ goo; and g++ = g The dif- 
ference between go± and goo is tiny and not resolved in this 
plot. Lengths and densities are expressed in units of r c . 



neutral solvent species. For (a) and (b) respectively, 10 6 
and 10 5 MC configurations were required to reduce the 
errors to an acceptable level at large r (note that there 
are ten times as many particles for the latter state point). 
A box of size (20er) 3 is necessary to reach out to r « 9er. 
Each state point required nearly 3000 hours of CPU time, 
in marked contrast to the HNC solution which takes less 
than a second. This is the origin of the claim earlier that 
HNC can be up to ten million times faster than MC. 

MC simulation of the thermodynamics is much less 
demanding. Fig. [2] shows the excess pressure (p cx = 
p — pk-^T) and internal energy for the URPM along 
isotherms at Ib/v = 1 and 10. There is excellent quan- 
titative agreement between HNC and MC. For the most 
part these simulations were carried out in a box of size 
(10cr) 3 . For Ib/& = 1 though we did check for finite size 
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FIG. 4. Total correlation function h++(r) from HNC for the 
URPM at l B = 1 and p z = 0.01(1)5 (top to bottom). Curves 
have been displaced for clarity. The Kirkwood transition from 
pure exponential decay (p z < 0.03) to damped oscillatory 
(p z > 0.03) is clearly seen. Lengths and densities are ex- 
pressed in units of a. 
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FIG. 5. The Kirkwood line for the URPM. The solid line 
with circles is from HNC. The dashed line is the RPA, from 
Eq. ( |12[ |. If solvent particles and short range repulsions are 
added, this map is practically unchanged. 



effects at the lowest investigated density by increasing 
the box size to (15c) 3 and (20ct) 3 . The results are shown 
in Fig. [2]as the multiple data points at p 2 cj 3 = 0.01, and 
indicate that finite size effects are small. 

The excess pressure and the internal energy (divided 
by three) are both expected to trend to the Debye-Hiickel 
limiting law, p cx = (U)/3V — — fc 3 3 /(247r), as the density 
decreases. As can be seen, the approach is rather slow. 
Note that the relation p cx — (U) /3V is a consequence of 
Clausius' virial theorem applied to point particles inter- 
acting with the Coulomb potential [3D]. 

In the presence of a neutral solvent the attainable 
MC accuracy is much diminished, largely because of the 
need to equilibrate the solvent particles. Fig. [3j shows 
an example of pair distribution functions for a solvated 
model at a typical state point. Again there is excellent 
agreement between HNC and MC (box size (10r c ) 3 ). In 
this plot note that symmetry enforces go+ = go- an d 

g ++ = g The approximate symmetry go± ~ goo is 

only very weakly broken in HNC (and not at all in the 
RPA) since it is exactly true that U ± = U o- 

To summarise the key result of this section: HNC ac- 
curately reproduces MC for the parameter ranges of in- 
terest. Thus we conclude that we can place enough confi- 
dence in HNC to use it as a tool to explore the properties 
of the model. 



B. Screening properties from HNC 

We now use the HNC to calculate the screening length 
from the asymptotic behaviour of the computed total cor- 
relation functions, focussing at first on the URPM. Fig. [4] 
shows the typical behaviour of h++ (r) through the Kirk- 
wood transition as p z varies at fixed Ib ■ Note that the de- 
cay rate of the total correlation function at first increases 
with increasing density, until one crosses the Kirkwood 
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FIG. 6. (color online) The screening length for the URPM, 
comparing the value extracted by fitting the asymptotic tails 
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FIG. 7. (color online) Comparison between a fully solvated 
model (solid line, black) and the URPM equivalent (dashed 
line, blue). Both are calculated using HNC. The dotted line 
is the RPA prediction from Eq. (111. Lengths and densities 
are expressed in units of r c . 
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FIG. 8. (color online) HNC results for a 1:2 electrolyte. 
Lengths and densities are expressed in units of r c . 



line, after which the decay rate remains roughly constant 
but the period of the oscillations decreases. This be- 
haviour is similar to the RPA, and presumably reflects 
the pole structure as discussed in section |III C| 

One can 'zero in' on the Kirkwood line transition by 
systematically narrowing the range of densities which are 
plotted. In this case one finds the transition is located 
at p z er 3 sa (30.0 ± 0.5) x 10~ 3 . By proceeding in this 
way, the entire Kirkwood line can be mapped out in the 
(p z cr 3 , Zb/c) plane. This is shown in Fig.fsl where HNC is 



compared to the RPA Kirkwood line from Eq. (12). We 
see that for Zb/ct < 5 the RPA is practically indistinguish- 
able from HNC, and even at — 10 the difference is 
still less than 40%. Above this value of Zb/ct HNC ceases 
to converge to a solution. 

On the low density side of the Kirkwood line the HNC 
screening length can be found by fitting the tail of to- 
tal correlation function to the expected asymptotic be- 
haviour in Eq. Q. Results along isotherms at three val- 
ues of Zb/c ar e shown in Fig. |6j Crucially, we see that 
for Zb/c < 5, the RPA screening length from Eq. ( 11 ) is 
in error by less than 10% compared to HNC. 

All the results presented so far have been for the 
URPM in the absence of a neutral solvent species. Re- 
markably, we have found that very little changes if short 
range repulsions are added (A — 25) and a solvent is 
included (pr 3 = 3). For example the Kirkwood line in 
Fig. [5] is practically unchanged and we have found the 
same to be true for the screening length itself. We give a 
single example here. Fig. [7] shows the asymptotic decay 
of h ++ (r) for the indicated state point for a fully solvated 
model, compared to the equivalent URPM at Zb/c = 2 
and p z a 3 = 0.0125. A line indicating the RPA decay 



from Eq. (11 1 is also included. We see that the presence 
of a solvent and short range repulsions confers some liq- 
uid structure at short distances but the asymptotic decay 
is practically unchanged, and agrees well with the RPA. 

Lastly we turn to the more complicated case of an 
asymmetric electrolyte. Fig. [8] shows the total correla- 
tion functions for a 1:2 electrolyte calculated using HNC. 
As can be seen the asymmetry splits apart the three ionic 
correlation functions, but nevertheless they share a com- 
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FIG. 9. (color online) Ratio between RPA screening length 
and Debye length for a 1:1 electrolyte, as a function of concen- 
tration, for three choices of a. The lower (upper) horizontal 
axis shows the concentration in physical (simulation) units. 
Each curve terminates when the model system crosses the 
Kirkwood line. The dashed line is at Arpa/Ad = 0.95. 
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TABLE I. Sample calculation for a 0.1 M 1:1 electrolyte. 



mon decay length, AhncAc ~ 0.93. This can be com- 
pared to Arpa/^c « 1.02, calculated from Eq. (11 ) using 
Pz r c = 0-06 (i. e. p z = ( o + +4/7_). The difference between 



HNC and RPA is less than 10%, as might be expected 
from Fig.|6]since the 1:2 case is intermediate between the 
1:1 case (Zb/ct = 2) and the 2:2 case (Zb/ct = 8). 

The main conclusions from this section are: first the 
solvent has practically no effect on the screening proper- 
ties so that Fig. [5] can be used as a quasi- universal quide, 
and second for many applications, such as to 1:1 elec- 
trolytes, the RPA suffices. 



C. Worked example 

Let us work through the example given at the end of 
section |IIIB[ for which the RPA solution is applicable. 
The calculations are shown in the numbered rows in Ta- 
bic [I] We start from the standard DPD mapping with 
r c = 0.645 nm and Zb = 0.7nm. This gives Zb/t" c = 1.09 
(row 1). If the molar concentration of a 1:1 electrolyte is 
c s , then p z = 2 x 10 3 c s iV^i (row 2; the factor two accounts 



for both species of ion). The Debye length (row 3) fol- 
lows from Eq. ([5]), here in the form Ad = (4-kIbPz) 1 ^ 2 ■ 
Alternatively one can use the well known expression 
Ad = O.Slnm/y^cJ from the colloidal literature [3T]. We 
choose a value for a (row 4) and calculate the left hand 
side of the inequality in Eq. ( [l2| (row 6). For the choice 
a = r c /2, Eq. (12) is satisfied and we are on the low 
density side of the Kirkwood line. We can then use the 
value of the Lambert function (row 7) to calculate Arpa 
(row 8) from Eq. ( 11 ). Given that Ib/c < 5 (row 5), this 



should be a good estimate of the true screening length. 
As claimed in section [IlIB[ the answer deviates from the 
Debye length by only 6% (row 9). For the choice a = r c 
though, Eq. (12 1 is violated (row 6), and we are on the 
high density side of the Kirkwood line. This is also in- 
dicated by the fact that the Lambert function (row 7) 
evaluates to a complex number. 



D. The choice of a 

As we have seen, a mapping to a physical system fixes 
Ib and r c , but the choice of a remains unresolved. Our 
study so far reveals this choice is a balance of conflicting 
requirements. On the one hand we would like to increase 
a as much as possible, mainly because this reduces the 
cost of computing the electrostatic interactions in a simu- 
lation. On the other hand if a is too large we run the risk 
of deviating strongly from the expected screening proper- 
ties of the physical system, and may ultimately cross the 
Kirkwood line in Fig. [5] Such behaviour is almost cer- 
tain to be artefactual since the chances of coinciding with 
similar behaviour in the physical system seem remote. 

For example for a 1:1 aqueous electrolyte we can plot 
the ratio Arpa/Ad as a function of salt concentration 
c s , using the method just described in section |IV C| 
Fig. [9] shows just such a plot, for three choices of a/r c . 
Inspection suggests as sensible compromise might be 
a/r c = 0.5, since this restricts significant deviations from 
the Debye length (i. e. more than 10%) to c s > 0.15 M, 
where in any case the Debye length is starting to become 
comparable to r c . 



V. DISCUSSION 

Let us close with some more remarks about implemen- 
tation, and indicate avenues for future work. First, let 
us dispose of an elementary point. The simulations de- 
scribed here have been performed using MC, rather than 
DPD. The reason for this is that we are interested in 
equilibrium properties, and MC is free from issues such 
as the choice of integration algorithm and time step [3j. 
Nevertheless the Ewald method can easily be applied to a 
dynamical simulation, by calculating the forces that arise 
from the potential energy in Eq. (13). 



The usual Ewald implementation for point charges in- 



troduces a 'splitting parameter' so that part of the inter- 
action is calculated in real space and part in reciprocal 
space [28]. Here we have 'physical- ised' the splitting pa- 
rameter by linking it to the Gaussian charge size a, so 
that we can discard the real-space interaction. This may 
not always be the best choice, since one cannot then op- 
timise the splitting parameter to match the simulation 
box size [jQ . However Coslovich et al. found that there is 
practically no benefit in divorcing the splitting parameter 
from the Gaussian charge size, at least for the URPM for 
their parameter ranges. Nevertheless it is worth bearing 
in mind this possibility, particularly if a is much smaller 
than the simulation box size. 

Aside from standard Ewald, any existing molecular dy- 
namics (MD) method could be used in principle to calcu- 
late electrostatic interactions in DPD. Most notable are 
the P3M (particle-particle-particle-mesh) methods, such 
as that introduced by Groot [11], and hybrids such as 
smooth particle mesh Ewald [32]. Some of these meth- 
ods are highly parallelisable [33], or highly efficient in 
other ways [32 . These MD methods are typically devel- 
oped for point charges, but the application to smeared 
charges should involve a straightforward extension to the 
underlying algorithms. 

As mentioned in the introduction, there is no consen- 
sus on the best form of charge smearing (linear, expo- 
nential, Gaussian, etc.), nor, perhaps, does there need to 
be. All smearing methods generate pair potentials which 
may differ in the short range part, but share a common 
Ib/t dependence for large r. This raises the question of 
whether the methods can be mapped on to one another, 
vis a vis the screening properties. Related to this is our 
observation that short range repulsions have practically 
no effect on the screening properties. Whilst this is a 
great bonus for applications, it cannot hold generally, for 
it would imply that there should be negligible effect of 
the choice of smearing. But we know this is manifestly 
untrue: two Gaussian charge models with a' 7^ a do not 
have the same screening properties. More generally, in 
any smeared charge model, another length scale must be 
present to non-dimensionalisc Ib- The question of how 
to determine this length scale remains unsolved. Our 
present results suggest that a systematic use of HNC 
could give an answer, thus providing a 'Rosetta Stone' 
tying together the existing treatments of DPD electro- 
statics. This is the subject of ongoing investigations. 

Separate from this, a long term goal is to incorporate 
specific ion effects into the model, such as the Hofmeis- 
ter series [33]. As mentioned in section [TlJ we have 
here focussed on a constant repulsion amplitude matrix 
Aij = A. Obviously there is scope to go beyond this, 
using HNC to calculate both the structural and thermo- 
dynamic consequences of unequal repulsion amplitudes. 
The hope is that a suitable choice of can be found, 
which will systematically and transferrably capture spe- 
cific ion effects. It is encouraging to note in this respect 
that a similar programme has been pursued with some 
success recently, for MD [M] [35] . 
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